loo
libraries
library(tidyverse)
library(brms)
library(bayesplot)
library(here)
library(glue)
library(scales)
theme_set(theme_light())
source(glue("{params$common_dir_str}/brms_model.R"))
source(glue("{params$model_dir_str}/model_prior.R"))data
obs_only <-
read_csv(glue("{params$model_dir_str}/data/stimulation_obvs.csv")) %>%
mutate(subj = as_factor(subj),
obs_degree = error,
error = obs_degree * (pi/180))## Parsed with column specification:
## cols(
## subj = col_double(),
## subj_index = col_double(),
## stimulation = col_double(),
## error = col_double()
## )full model
print(bprior_full)## prior class coef group resp dpar nlpar bound
## 1 normal(3.8, 0.5) b intercept circSD
## 2 normal(0, 0.5) b stimulation circSD
## 3 normal(0, 0.5) sd Intercept subj circSD
## 4 normal(0, 0.5) sd stimulation subj circSD
## 5 normal(0, 1.5) b intercept theta
## 6 normal(0, 1.5) b stimulation theta
## 7 normal(0, 1) sd Intercept subj theta
## 8 normal(0, 1) sd stimulation subj thetaiter = 4000
warmup = 2000
cores = 4
chains = 4
n_post_samples = (iter - warmup) * chains
fit_full <- brm(bform_full, obs_only, prior = bprior_full,
family = vm_uniform_mix, stanvars = stanvars,
warmup = warmup, iter = iter, cores = cores, chains = chains,
control = list(adapt_delta = 0.99), inits = 0,
file = glue("{params$save_dir_str}/fit_full"))fit check
divergences
model_fit <- fit_full
#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)
np %>%
filter(Parameter == "divergent__") %>%
summarise(n_div = sum(Value))## n_div
## 1 0rhat
mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
neff ratio
mcmc_neff(neff_rat) + yaxis_text(hjust = 1)
trace plots
mcmc_trace(as.array(model_fit$fit))
reduced model: no pre to post change on circSD parameter
print(bprior_DcircSD_null)## prior class coef group resp dpar nlpar bound
## 1 normal(3.8, 0.5) b intercept circSD
## 2 normal(0, 0.5) sd Intercept subj circSD
## 3 normal(0, 0.5) sd stimulation subj circSD
## 4 normal(0, 1.5) b intercept theta
## 5 normal(0, 1.5) b stimulation theta
## 6 normal(0, 1) sd Intercept subj theta
## 7 normal(0, 1) sd stimulation subj thetafit_DcircSD_null <- brm(bform_DcircSD_null, obs_only, prior = bprior_DcircSD_null,
family = vm_uniform_mix, stanvars = stanvars,
warmup = warmup, iter = iter, cores = cores, chains = chains,
control = list(adapt_delta = 0.99), inits = 0,
file = glue("{params$save_dir_str}/fit_DcircSD_null"))fit check
divergences
model_fit <- fit_DcircSD_null
#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)
np %>%
filter(Parameter == "divergent__") %>%
summarise(n_div = sum(Value))## n_div
## 1 0rhat
mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
neff ratio
mcmc_neff(neff_rat) + yaxis_text(hjust = 1)
trace plots
mcmc_trace(as.array(model_fit$fit))
reduced model: no pre to post change on pMem parameter
print(bprior_DpMem_null)## prior class coef group resp dpar nlpar bound
## 1 normal(3.8, 0.5) b intercept circSD
## 2 normal(0, 0.5) b stimulation circSD
## 3 normal(0, 0.5) sd Intercept subj circSD
## 4 normal(0, 0.5) sd stimulation subj circSD
## 5 normal(0, 1.5) b intercept theta
## 6 normal(0, 1) sd Intercept subj theta
## 7 normal(0, 1) sd stimulation subj thetafit_DpMem_null <- brm(bform_DpMem_null, obs_only, prior = bprior_DpMem_null,
family = vm_uniform_mix, stanvars = stanvars,
warmup = warmup, iter = iter, cores = cores, chains = chains,
control = list(adapt_delta = 0.99), inits = 0,
file = glue("{params$save_dir_str}/fit_DpMem_null"))fit check
divergences
model_fit <- fit_DpMem_null
#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)
np %>%
filter(Parameter == "divergent__") %>%
summarise(n_div = sum(Value))## n_div
## 1 0rhat
mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
neff ratio
mcmc_neff(neff_rat) + yaxis_text(hjust = 1)
trace plots
mcmc_trace(as.array(model_fit$fit))
reduced model: no pre to post change on pMem or circSD parameter
bprior_DcircSD_DpMem_null## prior class coef group resp dpar nlpar bound
## 1 normal(3.8, 0.5) b intercept circSD
## 2 normal(0, 0.5) sd Intercept subj circSD
## 3 normal(0, 0.5) sd stimulation subj circSD
## 4 normal(0, 1.5) b intercept theta
## 5 normal(0, 1) sd Intercept subj theta
## 6 normal(0, 1) sd stimulation subj thetafit_DcircSD_DpMem_null <-
brm(bform_DcircSD_DpMem_null, obs_only, prior = bprior_DcircSD_DpMem_null,
family = vm_uniform_mix, stanvars = stanvars,
warmup = warmup, iter = iter, cores = cores, chains = chains,
control = list(adapt_delta = 0.99), inits = 0,
file = glue("{params$save_dir_str}/fit_DcircSD_DpMem_null"))fit check
divergences
model_fit <- fit_DcircSD_DpMem_null
#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)
np %>%
filter(Parameter == "divergent__") %>%
summarise(n_div = sum(Value))## n_div
## 1 0rhat
mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
neff ratio
mcmc_neff(neff_rat) + yaxis_text(hjust = 1)
trace plots
mcmc_trace(as.array(model_fit$fit))
loo compare
#need to do this only once?
expose_functions(fit_full, vectorize = TRUE)
fit_full <- add_waic(fit_full, file = glue("{params$save_dir_str}/fit_full"))
fit_DcircSD_null <- add_waic(fit_DcircSD_null, file = glue("{params$save_dir_str}/fit_DcircSD_null"))
fit_DpMem_null <- add_waic(fit_DpMem_null, file = glue("{params$save_dir_str}/fit_DpMem_null"))
fit_DcircSD_DpMem_null <- add_waic(fit_DcircSD_DpMem_null, file = glue("{params$save_dir_str}/fit_DcircSD_DpMem_null"))
print("fit_full")## [1] "fit_full"fit_full$waic##
## Computed from 8000 by 504 log-likelihood matrix
##
## Estimate SE
## elpd_waic -804.9 15.9
## p_waic 7.2 0.2
## waic 1609.8 31.8print("fit_DcircSD_null")## [1] "fit_DcircSD_null"fit_DcircSD_null$waic##
## Computed from 8000 by 504 log-likelihood matrix
##
## Estimate SE
## elpd_waic -805.3 15.9
## p_waic 7.2 0.2
## waic 1610.5 31.7print("fit_DpMem_null")## [1] "fit_DpMem_null"fit_DpMem_null$waic##
## Computed from 8000 by 504 log-likelihood matrix
##
## Estimate SE
## elpd_waic -804.1 15.9
## p_waic 6.5 0.2
## waic 1608.3 31.8print("fit_DcircSD_DpMem_null")## [1] "fit_DcircSD_DpMem_null"fit_DcircSD_DpMem_null$waic##
## Computed from 8000 by 504 log-likelihood matrix
##
## Estimate SE
## elpd_waic -804.6 15.9
## p_waic 6.7 0.2
## waic 1609.3 31.7loo_compare(fit_full, fit_DcircSD_null, fit_DpMem_null, fit_DcircSD_DpMem_null, criterion = "waic")## elpd_diff se_diff
## fit_DpMem_null 0.0 0.0
## fit_DcircSD_DpMem_null -0.5 0.6
## fit_full -0.8 0.2
## fit_DcircSD_null -1.1 0.7*Notes
All models were fit well.
WAIC diffs were not very large and had non-neglible SE. Also inconclusive.
waic diff ordering:
fit_DpMem_null
fit_DcircSD_DpMem_null
fit_full
fit_DcircSD_null
fit_DpMem_null predicted best since it omits the difference parameter on pMem, for which there seems little evidence of an effect.
fit_DcircSD_null fit worst since it omits the parameter that seems to be actually showing a non-zero value.
fit_full and fit_DcircSD_DpMem_null are in the middle.